Towards a quantitative phase-field model of two-phase solidification 
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We construct a diffuse-interface model of two-phase solidification that quantitatively reproduces 
the classic free boundary problem on solid-liquid interfaces in the thin-interface limit. Convergence 
tests and comparisons with boundary integral simulations of eutectic growth show good accuracy 
for steady-state lamellae, but the results for limit cycles depend on the interface thickness through 
the trijunction behavior. This raises the fundamental issue of diffuse multiple-junction dynamics. 
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Complex microstructures that arise during alloy solid- 
ification are a classical example of pattern formation [Ij 
and influence the mechanical properties of the finished 
material 0. A long-standing challenge is to understand 
the pattern selection starting from the basic ingredients: 
bulk transport, solute and heat rejection on the solidifi- 
cation front, and the front's local response. Simple as it 
may seem, this free boundary problem (FBP) accurately 
describes many experimental features, but has few ana- 
lytic solutions, so that numerical modeling is mandatory. 

The phase-field method |3| has become the method of 
choice for simulating solidification fronts 4], and more 
generally for tackling FBPs and interfacial pattern for- 
mation phenomena, e.g. in materials science and 
fluid flow 0. Its main advantage (essential in three di- 
mensions) is that it circumvents front tracking by using 
phase fields to locate the fronts. These fields interpo- 
late between different constant values in each bulk phase 
through interfacial regions of thickness W . The model 
is then required to reproduce the FBP in the sharp- 
interface limit, in which the extra length scale W van- 
ishes. 

In practice, simulations have to resolve the variation of 
the phase fields through the interfaces, so that W must 
stay finite. Their results generally depend on the ratio 
W/£, where £ is a relevant length scale of the FBP. Ex- 
plicit corrections to the original FBP to first order in W/l 
have been calculated by a so-called thin- interface analy- 
sis in a few cases, and some canceled out [fj, LD, la, l2| • A 
complete cancellation, achieved for single-phase solidifi- 
cation 0, EJ , means that results become independent of 
W/t for some finite value of W. The correct FBP is then 
reproduced already at that value, much larger than the 
thickness of real interfaces, enabling quantitative contact 
in three dimensions between simulations, theory, and ex- 
periments in reasonable simulation times |lfj|. 

Here, we extend these advances to two-phase solidifica- 
tion, which already includes the most widespread solid- 
ification microstructures after dendrites: eutectic com- 
posites. They consist of alternate lamellae of two solids 
(a and (3) or of rods of one solid embedded in the other, 
growing from a melt L near a eutectic point, where all 
three phases coexist at equilibrium. The interplay be- 



tween capillarity and diffusive bulk transport between ad- 
jacent solid phases can give rise to more complex patterns 
and nonlinear phenomena such as bifurcations, limit cy- 
cles, solitary waves, and spatiotemporal chaos 

A two-phase solidification front consists of (i) solid- 
liquid interfaces and (ii) trijunction points where all three 
phases meet. Our strategy is to construct a phase-field 
model that allows us to analyze the thin-interface behav- 
ior of (i) separately from (ii). We quantitatively repro- 
duce the correct FBP on (i); (ii) satisfy Young's law at 
equilibrium. We test convergence in W/ £ for lamellar eu- 
tectic growth at experimentally relevant parameters, and 
compare our results to boundary integral (BI) jl2 S| simu- 
lations and other phase-field models. For steady states, 
we achieve good agreement with the BI and a drastically 
improved, fast convergence compared to previous mod- 
els. In contrast, convergence is slow for limit cycles, due 
to a trijunction behavior affecting the overall dynamics. 

We use one phase field pi to indicate presence (pi = 1) 
or absence (p, = 0) of each phase i = a, f3, L in the spirit 
of volume fractions , which requires 



Pa +P/3 +Pl = 1. 
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The phase fields evolve in time to minimize a free energy 
functional T of p = (j>ajP/3jPh), the solute concentration, 
and temperature, 
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Pq+P/3+Pl = 1 



where r(p) is a phase-dependent relaxation time. This 
classical problem of minimizing a functional subject to a 
constraint is treated by the method of Lagrange multi- 
pliers; {8T/5pi)\ Pct+P(i+Plj=1 = 5F/8pi-{\/$) Ysj ST/Spj 
for three phases, where the functional derivatives on the 
r.h.s. are now taken as if all pi were independent. 

To distinguish between phases, earlier phase-field mod- 
els of two-phase solidification used either the usual solid- 
liquid phase field and the local concentration 01 or in- 
troduced a second, a-f3 phase field 0|. Across a solid- 
liquid interface, both fields must vary, so that their dy- 
namics are coupled, which complicates a thin-interface 
analysis. The same is true for a generic choice of T in 
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Eq. (J2J. However, if on an i—j interface we can assure 
that the third phase field pk is exactly zero, Pi or pj can 
be eliminated using Eq. JIJ, so that the interface can be 
described in terms of a single independent variable. This 
was rec ently achieved using a free energy with cusp-like 
minima |l6| , but no thin-interface analysis is available for 
that model. We also achieve absence of the third phase, 
but using a smooth free energy, by requiring pk = to be 
a stable solution for pk of Eqs. J2J) for each i-j interface: 
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Pa+P/3 +PL = l,Pfc=0 



The advantage is that the simplest choice for T yields a 
model that turns out to coincide with the quantitative 
model of Ref. 9] on those i—j interfaces. 

To construct our free energy, we split it into parts, 



F — I /grad + /TW + A/ c 
'V 



The first is a free energy penalty 
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for the gradients of the phase fields that provides the 
interface thickness W. The next is a triple-well potential 
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that generates the basic "landscape": one well per pure 
phase and "valleys" with double-well profiles along each 
Pk = cut, separated by a potential barrier on tri- 
junctions p a = pp = p\, = 1/3. The last part has a 
strength A (a constant that controls convergence) and 
couples the phase fields pi to the temperature T and the 
solute concentration C through c(C) = (C — Ce)/AC, 
where AC = Cp — C a , C a and Cp are the limits of the 
eutectic plateau, and (Ce,Ie) is the eutectic point, 



fc = Y,9i(p)lB i (T)-nA i (T)}, 



(7) 



where we have introduced the chemical-potential-likc 
variable /x = c — J2i Ai(T)hi, and g%{p) and hi(p) (given 
below) interpolate between for pi = and 1 for pi = 1 . 

The term f c drives the system out of equilibrium by 
unbalancing the pure phase free energies: Each well i is 
shifted by an amount Bi — \xAi. The equilibrium value 
A 4 = Mcq = {Bj ~ Bi)/(Aj — Ai) gives equal shifts and 
hence restores the balance between phases i and j; from 
the definition of /i, we obtain cf = Ai + /i*^ for the 



concentration in phase i coexisting with phase j. A eu- 
tectic phase diagram with constant concentration gaps 
and straight liquidus and solidus lines is generated by 
Ai = Ci = c(Ci) and Bi — a(T — Te)/ (m^AC), with rrii 
the (signed) liquidus slopes, i = a, (3. Non-constant con- 
centration gaps and peritectic phase diagrams can also 
be treated. Without loss of generality, A\, = £>l = 0. 

In order for [i — [i l J to keep the balance all across the 
i—j interface as pi goes from to 1, we require 

9i(pi,Pj,Q) = 1 - 9i(pj,Pi,0) Vi (8) 

Otherwise, several thin-interface corrections arise 0, . 
The simplest choice satisfying also Eq. (|3af) is gi = 



p?{15(l - Pi )[l + Pi - ipk-Pj?] +Pi(9p 2 - 5)}/4. 

The evolution of fi is obtained from its definition and 
mass conservation, dtc + V ■ J = 0, J = —Dp^Wfi + Jat- 



^ A at 



V • Jat, (9) 



where —Dpi,V^ is the usual diffusion current, with a dif- 
fusivity that varies from D in the liquid to in the solid 
(one-sided model), and Jat is an extension of the anti- 
trapping current introduced in that counterbalances 
spurious solute trapping, 

dpt 
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where hi — — Vpj/|Vpj| are unit vectors normal to i— L 
interfaces, and Ui • nL prevents solute exchange between 
the two solids. The model is not variational, because of 
the term Jat and because « ^ df c /dc, but enables us to 
use hi = Pi, which allows for a coarser discretization • 

Our model [Eqs. @ and ©] has stable interface solu- 
tions connecting two coexisting phases i and j: fi — fA^, 
Pi = 1 - pj = {1 ± tanh[r/(W^v / 2)]}/2 (with r the dis- 
tance to the interface), pk = 0. Since these solutions are 
identical for all i-j pairs, so are the i-j surface tensions. 
Unequal surface tensions can be obtained by adding new 
terms in Eq. J1J that shift the i-j free energy barriers. 

Remarkably, on solid-liquid (i—L) interfaces, assuming 
a weak dependence of the Ai, Bi on T, and r(p) = Tj, the 
change of variables </>,; = pi — Pl, u = (fj£ — fi)/Ai maps 
Eqs. J2J an< 4 © to the quantitative model with constant 
concentration gap in |9j, up to numerical prefactors. The 
thin-interface limit can hence be deduced by inspection 
and yields the classic FBP on i—L interfaces, 
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where Eq. I|lla|) holds in the liquid and the others are 
boundary conditions on the interface that has normal 



3 



velocity v n and curvature re; the minus (plus) refers to i = 
a (P) , and the capillary lengths di and kinetic coefficients 
Pi read in terms of our model parameters 
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with ai — v2/3 and a 2 = 1.175. The constant A oc W/di 
in Eqs. 10} , l |12(l and (|13|l controls the convergence to the 
original FBP. Any set of Pi can be treated with suitable 
Ti. We consider here p a — Pp = 0, which is achieved with 
n = a 2 AjXW 2 /D. The different n for A a ^ (e.g. 
different concentration gaps) are interpolated by r(p) = 
f+(l/2)(T a -T )(p a -p l3 )/(p a +p ), r(p a +p p = 0) =f, 
with f = (r Q + Tf})/2. 

We test our model in directional solidification with 
T = T E + G(z - Ft) , where G > is the thermal gradient 
and V > 0, the pulling speed, both directed along the z 
axis. Half a eutectic lamellae pair of total width A is sim- 
ulated in two dimensions [x and z) with no-flux boundary 
conditions in the midline of each lamella, using a finite- 
difference Euler scheme with a grid spacing Ax = 0.8W 
(coarser far into the liquid to improve efficiency). We 
adopt l D /d = 51200 and l T /l D = 4, where l D = D/V 
is the diffusion length, l T = \rrii\Ac/G are the thermal 
lengths, andJ= (d a +dp)/2, It = (l T +l T )/2. These cor- 
respond to typical experimental values G ~ lOOK/cm, 
V « lfim/s for CBr4-C2Clg, an organic eutectic for 
which accurate experimental data exist 0|- We use 
m a = —mp, c a — ~cp (a symmetric phase diagram) 
or mpl"m a — —2, —cp/c a — d a /dp — 2.5 (one close to 
CBr4-C2Cle). In both cases fi(z — > +oo) = (eutectic 
composition). We test convergence to the thin- interface 
limit with decreasing W by conversely increasing X/W 
while keeping all the ratios above and A/A m i n fixed, where 
Amin oc \J dlfj is the minimal undercooling spacing |l7j |. 
This is achieved by varying the constant A in Eq. (|12J) . 

Figure ^ shows the solid-liquid interfaces of a steady- 
state lamellae pair calculated by different phase-field 
models and the boundary integral method (BI) for 
A « A m i n . For the symmetric phase diagram [Fig. H^a)], 
our model (thin solid lines) agrees well with the BI (thick 
solid line). Moreover, the curves at X/W = 64, 92 and 
128 are indistinguishable. This means that the results 
are independent of X/W for X/W > 64, the signature of 
a quantitative model. In contrast, if we remove the anti- 
trapping current in our model, Jat = 0, which leads to 
solute trapping and finite interface kinetics, the results 
depend on X/W for all the range from 32 (bottom dashed 
line) to 128 (top one). The convergence of models not 
backed by a thin-interface analysis can even be slower, 
as shown by the dotted curves for a qualitative version 
of our model with hi = gi violating Eq. © and Jat = 
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FIG. 1: Steady-state lamellae pair profiles (dimensionless 
undercooling vs. x/X) for different models. Four curves at 
X/W =32, 64, 96 and 128 shown per model; curves closer to 
the boundary integral: larger X/W. [X/W = 64-128 collapse 
for the present model with antitrapping current in (a)]. Phase 
diagram used: (a) symmetric; (b) close to CBr4-C2Cl6. See 
parameters in the text. Inset: Averaged undercooling in (b) 
vs. X/W , compared to that without antitrapping current. 



|l8|: in this situation, several thin-interface corrections 
to the FBP occur simultaneously 0, . 

Results are similar for the phase diagram close to 
CBr4-C2Cle [Fig. ^b)]. The convergence is somewhat 
slower, since one of the lamellae is thinner and needs to 
be properly resolved. Some small deviation from the BI 
persists, probably due to the trijunction behavior (see 
below). In the inset, we plot the average undercooling 
vs. X/W. This is a less stringent test, as shown by the 
fact that results for our model are converged already for 
X/W — 32. However, those for the model with Jat = 
still depend on X/W at X/W = 128, which illustrates 
how all corrections need to be canceled before quantita- 
tive results can be achieved. 

Next, wg increase A to ~ 2 . 2 , close above the 
threshold A ps 2A min ^1 f° r the bifurcation from steady 
lamellae to oscillatory limit cycles, a situation in which 
the oscillation amplitude is very sensitive to all param- 
eters. Indeed, for the symmetric phase diagram and 
X/W = 64, the qualitative model of Ref. [18| still yields 
lamellae, whereas the present model correctly produces 
cycles, which are shown in Fig. [3(a) . However, the 
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FIG. 2: Limit cycles, (a) Superimposed snapshots of the 
interfaces at constant time intervals for X/W = 64. Thicker 
lines: a-/3 interfaces, (b) Amplitude of the trijunction oscil- 
lation in units of A vs. X/W. The line is a fit that yields 
A{X/W -> do) /A = 0.142. (c) Blowup of 6.4W x 6 AW.' Solid 
lines: trijunction passage; dashed line: later a-/3 interface. 



yield a substantial efficiency gain for small curvatures 
of trijunction trajectories, which makes it a promising 
tool for three-dimensional simulations. Second, the free 
boundary problem to converge to should also be recon- 
sidered. It was shown elsewhere that Young's condition 
on the angles between interfaces is violated out of equilib- 
rium for kinetically limited growth [19|; here, the global 
trijunction rotation was found to be fairly independent 
of the interface thickness, so that it might persist for real 
nanometric interfaces. These effects should be further 
investigated, possibly by atomistic simulations. 
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amplitude of the trijunction oscillation A/X, defined as 
its maximal displacement in x/X, strongly depends on 
X/W, as shown in Fig. EJb). An extrapolation yields 
A(X/W -> oo) = 0.142A, not far from the BI result 
A = 0.139A, but the results are still not converged for 
X/W = 192, in strong contrast to the steady-state behav- 
ior. This suggests that some correction(s) to the FBP in 
W/X remain in our model. Since solid-liquid interfaces 
are controlled, we turn to the the trijunctions. 

The solid (dashed) lines in Fig.|2fc) show a first (later) 
snapshot of the interfaces close to a turning point of the 
trijunction trajectory. In the later one the trijunction has 
moved away and only the a- (3 interface remains, which 
has slightly moved sideways. In the one-sided FBP, (i) 
the ot-(3 interface cannot move, so it is the trace left by 
the trijunction, and (ii) its direction close to the trijunc- 
tion approaches that of the trijunction velocity. In a 
diffuse-interface model, the diffusivity behind the trijunc- 
tion point p a = pp — pl = 1/3 falls to zero on the scale 
of W, so that (i) and (ii) do not hold. We consistently 
observe the displacement to be a fraction of W fairly 
independent of X/W, and the whole trijunction to be 
slightly rotated with respect to its velocity, features also 
observed for the steady state in Fig. ^d. This effect ex- 
plains the remaining mismatch between phase-field and 
BI in Fig.^J) and the slow convergence of A/X here. 

We have presented a phase-field model of two-phase 
solidification that coincides with the best models to date 
0,0 on solid-liquid interfaces, whose dynamics are com- 
pletely controlled. This has allowed us to identify the role 
of diffuse trijunctions in the convergence of the results. 
Understanding their dynamics is both a fundamental is- 
sue and a prerequisite for a fully quantitative modeling of 
multiphase solidification: First, a thin-interface analysis 
of the trijunction region in the phase-field model is lack- 
ing. Even so, our model is expected to be precise and 
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